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Summary 


Tho  overall  objective  of  this  research  and  development  (phase  I  and  the  proposed 
phase  II)  contract  is  to  develop  a  computer  program  for  calculating  the  flow  about 
aerodynamic  bodies  in  relative  motion.  The  technical  objective  of  the  Phase  I  effort 
was  to  develop  and  demonstrate  a  flow  analysis  procedure  capable  of  simulating  the 
unsteady  three-dimensional  compressible  viscous  flow  associated  with  a  turbine 
rotor/stator  configuration.  The  technical  question  of  how  best  to  perform  coupling 
between  computational  zones  moving  relative  to  each  other  was  addressed,  since  it 
has  the  most  significant  influence  upon  the  feasibility  of  the  flow  analysis  procedure. 

The  Phase  I  work  plan  consisted  of  four  main  tasks.  The  first  was  to  develop  a 
general  intprzone  boundary  condition:  a  procedure  to  couple  adjacent  zones  of  mesh 
in  relative  motion.  This  was  done  using  trilinear  interpolation  across  the  “sliding” 
boundary  dividing  the  moving  zones.  The  second  task  was  to  incorporate  moving 
mesh  terms  and  the  new  interzone  boundary  conditions  into  an  existing  three- 
dimensional  multi-zone  Navier -Stokes  code.  This  was  done  using  a  “layered”  soft¬ 
ware  structure  designed  to  minimize  the  amount  of  problem  specific  programming 
required.  The  third  task  was  t,o  demonstrate  the  analysis  on  a  turbine  rotor/stator 
flow  field  and  the  last  task  was  to  document  the  results. 

The  results  of  the  Phase  I  Research  and  Development  elfort  are  encouraging. 
They  have  shown  the  benefits  of  developing  a.  layered  software  structure,  as  well 
as  demonstrated  the  feasibility  of  the  coupling  procedure  for  interzone  boundary 
conditions  between  zones  in  relative  motion.  The  technical  feasibility  of  completing 
the  proposed  flow  analysis  code  is  high.  The  computer  program  resulting  from  the 
follow  on  funding  will  have  potential  applications  throughout  the  aerospace  industry 
foi  analysis  of  aerodynamic  bodies  in  relative  motion. 
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Nomenclature 


E 

F 

kt 

k 

n 

P 

P 

Pri 

Pr , 

PSi+ 1/2 

tf. 

R 
S  • 

u 

w 

t 

T 

Ui 

P 

7 

Rj 

P 

Pi 

4> 


Specific  heat  at  constant  pressure 
Specific  heat  at  constant  volume 
Internal  energy  per  unit  mass  (e  =  c„T) 

Total  Energy  per  unit  Volume  E  =  p[e  +  1/2 (u2  +  v2  4-  u/2)] 

Fluxes  in  the  i- component  direction 

Laminar  coefficient  of  heat  conduction 

Coefficient  of  heat  conduction  (molecular  and  turbulent) 

Unit  outward  normal  vector  to  a  closed  surface,  S,  in  equation  1 
Static  Pressure 

Combined  flux  vector, P  =  Fp  -+-  F2j  +  P3fc,  in  equation  1 

Laminar  Prandtl  number,  Pr;  =  =  0.72 

Turbulent  Prandtl  number,  Prt  =  — ^  =  0.9 

Fluxes  through  surface  separating  cells  i,  j,  k  and  i  -+-  1 ,  j,  k 

Combined  molecular  and  turbulent  heat  fluxes 

Gas  Constant 

A  closed  surface  surrounding  volume,  V,  in  equation  1 
Conservative  Variables,  U  =  [p,  pu,  pv,  pw,  E]T 
Primative  Variables,  W  =  [u,u,u>,  p,  T]1 
Time 

Temperature 

ith  component  of  velocity 

1  CW.l^onuiL  Ui  puoiti,CMA 

Kronecker  delta  =  1  when  i  —  j  and  6 ,j  =  0  otherwise) 
Density 

Ratio  of  specific  heats 

Combined  viscous  and  turbulent  stress  tensor 
Combined  molecular  and  turbulent  viscosity  coefficient 
Laminar  viscosity  coefficient 
Parameter  in  interpolation  to  surface 
(< p  —  0;  First  order, ^  —  1,  Second  order) 

P aramctcr  in  interpolation  to  surface 

(k  =  —1;  Fully  upwind  differencing,  k  =  1;  Central  differencing) 


2  Introduction 


Computer  codes  that  numerically  simulate  the  flow  of  fluids  are  increasingly  be¬ 
ing  used  in  the  design  and  analysis  of  aircraft  components.  Although  full  three- 
dimensional  viscous  calculations  can  be  expensive,  they  also  can  provide  informa¬ 
tion  that  is  extremely  beneficial  Use  of  flow  analysis  codes  can  reduce  costs  of 
parametric  testing  and/or  increase  the  variety  and  extent  of  configurations  that  are 
considered  during  a  design  process.  Many  times  a  component  must  be  designed  for 
conditions  outside  any  existing  experimental  data  base  or  at  conditions  that  make 
parametric  testing  impractical.  For  such  situations  numerical  simulation  may  be 
the  principal  method  used  in  the  design  process. 

The  design  of  aircraft  components  has  evolved  to  a  sophisticated  level.  Many 
critical  components  such  as  the  flow  passages  in  turbomachinery,  nozzles,  and  inlets 
operate  at  very  high  efficiencies  (at  least  at  their  design  points).  Further  perfor¬ 
mance  improvements  will  be  difficult  and  expensive  to  obtain  using  'cut  and  try’ 
testing  methods  or  even  two-dimensional  flow  analysis  procedures.  It  appears  that 
three-dimensional  flow  analysis  capable  of  handling  complex  configurations  will 
eventually  provide  the  only  economical  method  for  designing  aircraft  components 
that  have  significantly  greater  performance. 

One  area  of  particular  interest  to  the  U.S.  Army  is  that  of  unsteady  flows  in 
which  one  or  more  aerodynamic,  bodies  move  relative  to  a  fixed  structure.  Examples 
of  this  flow  situation  axe:  turbomachinery  flow  passages,  helicopter  blade /fuselage 
interaction,  propeller /wing/nacelle  interactions,  dispensing  of  submunitions  from 
missiles,  missile  stage  separation  transient,  and  release/launch  of  munitions  from 
aircraft.  All  of  these  flow  situations  have  complex  geometries,  are  inherently  un¬ 
steady,  and  have  significant  viscous  effects.  The  experimental  data  of  Bring,  et 
al[8]  show  that  the  temporally  var  ying  pressure  measured  at  the  leading  edge  of  a 
turbine  rotor  can  be  as  much  as  72%  of  the  exit  dynamic  pressure  when  the  spac¬ 
ing  between  moving  rotor  and  fixed  stator  blades  is  small.  Obviously,  the  unsteady 
effects  in  turbom achinery  can  be  substantial.  There  is,  therefore,  a  need  for  an  engi¬ 
neering  analysis  tool  for  simulating  unsteady  viscous  flows  past  three-  dimensional 
aerodynamic  bodies  in  relative  motion. 

Numerous  calculations  of  cascade  flows  have  been  described  in  the  literature 
(for  example[9,l(J,5,7]).  Except  for  the  work  of  Rai[5,7]  all  of  the  known  work  has 
been  steady  flows  or  has  used  equation  sets  simpler  than  the  Navier-Stokes  equa¬ 
tions.  The  work  by  Rai  in[7]  suggests  that  the  flow  analysis  technology  is  available 
and  is  practical  for  full  three-dimensional  unsteady  turbomachinery  flow  calcula- 
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tions.  However,  all  of  the  existing  computer  codes  for  cascade  flows  are  ‘research 
codes’-codes  written  to  develop  algorithms.  They  are  not  engineering  tools.  An 
engineering  tool  needs  to  be  easy  to  use.  It  must  be  reliable  and  have  a  minimum 
number  of  input  variables  that  must  be  ‘tuned’  (e.g.  smoothing  terms)  to  achieve 
a  satisfactory  solution.  It  must  include  a  mesh  generator  for  at  least  one  specific 
class  of  configurations.  It  must  have  comprehensive  and  readable  documentation 
of  its  use,  limits  of  applicability,  and  theory.  It  must  be  thoroughly  validated.  In 
addition,  for  it  to  be  user  friendly,  it  must  check  the  user  input  data  for  errors  and 
then  give  informative  messages  back  to  the  user.  The  real  benefit  of  flow  analysis 
in  the  design  of  airciaft  components  will  only  be  realized  after  engineering  tools  of 
the  sort  described  above  become  a  nilable  to  design  engineers. 


This  contract  (Phase  I  and  the  proposed  Phase  II)  would  lead  to  an  engineering 
tool  for  the  numeneal  simulation  of  three-  dimensional  flows  around  aerodynamic 
bodies  moving  past  fixed  structures.  This  flow  analysis  procedure  would  divide  the 
computational  domain  into  one  or  more  zones  in  which  a  single  body-fitted  mesh 
can  easily  be  generated.  Mesh  zones  would  also  be  allowed  to  move  in  time-each 
with  a  unique  velocity  and  acceleration  to  provide  capability  for  moving  bodies.  The 
Phase  I  effort,  being  reported  here,  was  intended  to  determine  the  feasibility  of  the 
engineering  tool  discussed  above.  The  Phase  I  work  concentrated  a  single  problem — 
the  unsteady  rotor/stator  interaction  problem  solved  previously  by  Rai[5,7].  The 
work,  plan  consisted  of  the  following  four  tasks 


Develop  and  Incorporate  Interzone  Bou  :  ’.ary  Conditions  into  an  Existing 


nr  1 1  •  nTV  \T _ Pl..1 _ /'i  _  J  _ 

muiLi-ZiUiit;  ndviei-juu\c&  \juuc 


2.  Incorporate  Moving  Mesh  Terms  into  the  Navier-Stokes  Code 

3.  Demonstrate  the  Analysis  on  a  Turbine  Rotor/Stator  Flow  Field 

4.  Document  Results 


The  work  on  Task  1  is  described  in  detail  in  Section  3.2.2, 
described  in  Section  5. 


and  Task  3  is 


5 


3  Solution  Procedure 


The  pilot  code  solves  the  compressible  Navier-Stokes  equations  using  a  time-marching 
solution  procedure  and  a  multiple  zone  grid.  The  solution  procedure  is  discussed 
in  this  section  and  the  interzonc  boundary  conditions  are  discussed  in  Section  3. 
The  discussion  of  the  solution  procedure  is  broken  into  three  parts:  Section  2,1 
describes  the  mathematical  model,  Section  2.2  describes  the  solution  procedures 
used,  and  Section  2.3  describes  the  boundary  conditions  used.  Innovations  devel¬ 
oped  during  Phase  I  include  the  moving  mesh  terms  described  in  Section  2.2  and 
the  periodic  boundary  conditions  described  in  Section  2.3.  The  principle  innovation 
is  the  general  sliding  boundary  described  in  Section  3.2. 


3.1  Mathematical  Model 


The  governing  equations  that  are  solved  by  the  pilot  code  axe  the  time-dependent 
three-dimensional  compressible  Navier-Stokes  Equations.  These  equations  are  given 
below  in  integral1  form[l], 

£  guw  +  gpjtds-o  (1) 


where 

and 


P  =  F\i  -j-  F ij  +  F$k 


U 


F, 


'  P  N 
PH  l 
pu2 

PU'3 

\  E  / 

''  pui  > 

pu,ux  +  pSu  -  r,  1 
pu,u2  +  p62i  -  r. 2 
puiu3  +  pS3i  -  r,3 

\  (E  T  P)ui  -  ujru  +  / 


where  t  is  the  time,  V  is  an  arbitrary  volume,  5  is  the  surface  surrounding  V.  In  the 
above  equations  the  standard  summation  convention  (sum  over  repeated  indices)  is 
used. 

lThe  integral  form  of  the  equations  is  strongly  conservative. 


6 


The  following  auxiliary  equations  are  needed  to  close  the  system. 


E 


<U 


P 

e 


P 

P 

-k 


e  +  -ujuj 


1 
2 
du, 
dx 
8T 


duj 

lAdsj  +  dx, 


dx, 

p  (7  —  1)  e  =  pRT 
c„T 


2  duk 

r,Jdxk 


Here  p  and  k  are  the  sums  of  the  molecular  and  turbulent  values  of  the  viscosity 
and  heat  conduction  coefficients. 


The  molecular  viscosity  is  obtained  from  Sutherlands  law  for  air. 


pi  =  1.4519(10)~6 


^3/2 

r  + 110.0 


The  coefficient  of  heat  conduction  is  then  obtained  by  assuming  a  constant  Prandtl 
number,  Pr;  =  0.72. 


For  turbulent  flows  the  algebraic  model  of  Baldwin  and  Lomax [2]  is  used  to  cal¬ 
culate  pt.  The  turbulent  heat  conduction  coefficient  is  then  calculated  assuming  a 
turbulent  Prandtl  number  of  0.9. 

Physically  Equation  1  represents  a  very  simple  idea:  the  time  rate  of  change 
of  mass,  momentum,  and  energy  within  an  arbitrarily  chosen  volume,  V,  is  equal 
to  the  apparent  flux  of  these  quantities  inward  through  the  surface,  S,  surrounding 
the  volume.  The  finite  volume  method  consists  of  breaking  the  flow  field  up  into 
a  large  number  of  nearly  hexahedral  finite  volume  cells,  as  shown  in  Figure  1,  and 
applying  the  integral  equations  directly  to  each  volume.  The  finite  volume  method 
is  the  natural  form  for  a  conservative  solution  procedure  and  it  avoids  certain  metric 
singularities  which  require  special  treatment  in  a  finite  difference  method. 
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Figure  1:  Finite— Volume  Mesh 


3.2  Solution  Procedure 

The  pilot  code  is  a  zonal  code  using  body-fitted  nonorthogonal  grids  for  each  zone. 
Each  zonal  grid  is  structured  such  that  is  has  a  logical  i,j,  k  ordering.  Using  a 
large  number  of  zones  patched  together  a  body-fitted  grid  may  be  generated  for 
very  complex  geometries.  The  zonal  grids  are  generated  separately  and  read  into 
the  pilot  code  from  a  file.  The  user  specifies  the  zonal  connectively. 

The  Navier-Stokes  equations  are  solved  using  a  time-maxching  finite-volume 
method.  The  method  is  derived  by  applying  the  integral  equation,  Equation  1,  to 
a  finite-volume  given  by  indicies  i,j,  k  (see  Figure  2). 

jt  (UtihkVol,,hk)  =  -  [D,P-S  +  DJP-S  +  D,P  S)^k 


where  Ui,]tk  is  the  mean  value  of  U  in  cell  i,j,  k  and  DiP-S  for  example,  represents 
the  difference  of  the  fluxes  through  opposing  faces  of  the  cell  in  the  i-index  direction. 

The  time  derivative  is  approximated  using  forward  in  time  differencing: 


S 


Figure  2:  Finite- Volume  Cell 


Vol,, 

At 


=  -  [DiP-S  +  DjP-S  +  DkP-S )*\ . 


(2) 


where 


Wu*  =  psj  - 1^* 

The  fluxes  must  now  be  approximated  in  terms  of  the  Ui,j,k.  For  conciseness  consider 
only  the  i  +  1/2  surface. 

For  the  approximation  of  the  flux  through  a  surface  the  inviscid  and  diffusion 
terms  of  the  flux  vector  are  considered  separately. 


nmo  .  t~»  r*di  f  f 

r  •  o  -t-  r  ■  o 


These  terms  sue  then  evaluated  in  a  manner  consistent  with  the  predominant  nature 
of  the  equations  in  the  limit  as  iie  — ►  oo  (hyperbolic)  and  Re  — ►  0  (parabolic);  i.e., 
upwind  differencing  for  the  inviscid  terms  and  central  differencing  for  the  diffusion 
(viscous  stress  and  heat  flux)  terms. 
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3.2.1  Inviscid  Flux 


The  evaluation  of  the  inviscid  terms  is  based  on  the  flux  vector  splitting  method 
developed  by  Steger  and  Warming. [3]  This  method  splits  the  llux  vector  P'S  into 
vectors  containing  the  information  propogating  in  the  direction  of  S,f+1  and  in  the 
direction  opposite  to  S,f~.  The  flux  vector  for  the  i-f  1/2  surface  is  then  written 

The  notation  /T(I7~)  denotes  /+  evaluated  at  U~.  U~  and  U+  represent  upwind 
interpolations  to  the  ceil  interfaces 


U,ll/2  =  Ui.hk  ±  £  [(1  T  «)  Vi  +  (1  ±  k)  A,]  Uij,k  (3) 

where  YtU>,;,k  =  f —  t/,-_ i.j^,  is  a  backward  difference  operator  in  the  i-direction 
and  A,Uit}ik  —  is  a  forward  difference  operator  in  the  i-direction. 

The  parameter,  k,  determines  the  nature  and  accuracy  of  the  spacial  differencing: 
k  =  —  1  corresponds  to  a  second-order  fully  upwind  scheme,  k  =  +1  to  a  second- 
order  central  difference  scheme,  and  k  =  1/3  to  a  third-order  upwind  biased  scheme. 
For  first-order  accuracy  4>  is  set  to  zero  and  for  higher  order  accuracy  <p  is  set  to 
one. 


In  general,  some  form  of  flux  limiting  is  required  to  avoid  oscillations  at  shock 
waves.  The  limiting  can  be  implemented  by  locally  varying  the  difference  quotients 
in  the  interpolation,  Equation  3,  in  re’ation  with  local  gradients  of  the  solution. 
If  properly  implemented  a  limiter  will  yield  sharp  monotonic  shock  waves.  In  this 
investigation  a  simple  form  of  limiting,  based  on  the  second  difference  of  pressure, 
was  used.  In  particular 

4>ik\! 2  =  MAX  (1  -  CltmDDP,  0.0) 

where  DDP  is  the  magnitude  of  the  second  difference  of  pressure  normal  to  the 
surface  (in  this  case,  the  f-direction).  The  differencing  parameter,  «,  remains  un¬ 
changed.  This  limiting  reduces  the  form  of  differencing  toward  first-order  upwind 
differencing  in  regions  of  strong  pressure  gradients  (such  as  shock  waves).  The  spa¬ 
cial  order  of  accuracy  remains  second-order  everywhere  unless  is  actually  zero. 
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3.2.2  Moving  Mesh  Fluxes 


The  intregal  equation,  Equation  1,  may  be  applied  to  a  finite  volume  mesh  that  is 
moving.  The  limits  of  the  integrals  are  then  functions  of  of  time  and  their  eval¬ 
uation  is  more  complicated.  In  the  case  where  the  volumes  of  the  cells  remain 
constant  -rigid  body  motion  of  a  zone-the  descrete  conservation  equation,  Equa¬ 
tion  3.2,  retains  the  same  form  and  the  only  affect  of  moving  the  mesh  is  to  modify 
the  fluxes,  PS.  Ti  e  modified  fluxes  are 


(P-S)m  -  vnU\S | 


where  ( P-S)m  is  the  flux  function  without  moving  mesh,  vn  is  the  velocity  of  the 
surface  normal  to  itself,  and  U  is  some  mean  value  of  the  conservative  variables. 
During  Phase  I,  U  was  taken  as  the  average  of  the  conservative  variables  in  the  cells 
on  either  side  of  the  surface. 


3.2.3  Viscous  Flux 


The  contribution  to  the  fluxes  from  the  diffusion  terms  (viscous  stresses  and  heat 
conduction)  are  approximated  using  standard  central  differences.  Following  the  de¬ 
velopment  in  Reference^  1]  these  fluxes  can  be  written  in  terms  of  seco.  d  differences 


■\f  f  Ti  <o 

y k.  iuiv  UA  iitiuvi  *  U 


variables  in  the  three  generalized  coordinate  directions 


P-SliW  =  Fd(W(,W,n  (4) 

where  £,  t /,  and  £  correspond  to  the  directions  of  increasing  i,  j,  and  k  indices, 
respectively.  These  derivatives  are  approximated  as  follows 


Wf. 

Wv 

wc 


Wi+li]ik  -  WiJJt 

0-25  (TF,+l  j+1r  —  l^i+ip-U:  +  Wi,j+i,k  ~ 
0.25(Wi+1iM.+1  -  Wi+ljifc_1  +  -  Wij,k. i) 


(5) 


Using  equations  5  and  the  relationship  between  the  conservative  variables  and  the 
primative  variables  the  discrete  approximation  to  the  diffusion  fluxes,  Equation  4,  is 
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written  in  terms  of  the  conservative  variables  within  the  10  cells  nearest  the  i+1/2 
face.  The  thin  layer  form  of  the  equations  are  obtained  by  neglecting  Wt|  and 
in  Equations  4  and  5. 


3.3  Boundary  Conditions 

Five  types  of  boundary  conditions  are  needed  for  the  rotor-stator  interaction  prob¬ 
lem:  subsonic  inflow,  subsonic  outflow,  noslip  wall,  periodic,  and  interzone.  The 
first  four  boundary  conditions  will  be  discussed  here  while  discussion  of  the  last  wifi 
be  deferred  to  Section  4. 

The  implementation  of  ail  boundary  conditions  is  simplified  by  the  use  of 
boundary  cells.  Those  are  extra  cells  added  outside  the  solution  domain  solely 
for  the  purpose  of  satisfying  boundary  conditions.  When  boundary  cells  are  used 
the  solution  process  within  the  interior  of  the  solution  domain  may  proceed  with¬ 
out  any  special  differencing  near  the  boundaries.  The  current  solution  procedure 
reuquires  two  layers  of  boundary  cells. 


3.3.1  Subsonic  Inflow 

The  treatment  of  the  inflow  boundary  conditions  is  guided  by  the  theory  of  charac¬ 
teristics.  A  locally  one-dimensional  flow  has  five  characteristic  equations  with  slopes 
u'\  u'\  u\  u'  -}-  c,  and  u'  -  c.  If  the  flow  field  is  supersonic  then  all  five  characteristic 
equations  are  propogating  information  in  the  positive  x-direction.  In  this  case  all 
data  must  be  specified  at  the  inflow  boundary.  If  the  flow  is  subsonic:  at  the  inflow 
boundary,  then  one  of  the  characteristics,  the  u'  —  c  characteristic,  has  a  negative 
slope  and  it  propogates  information  from  the  interior  upstream  to  the  inflow  bound¬ 
ary.  In  this  case  only  four  items  may  be  specified  at  the  inflow  boundary  and  the 
fifth  item  must  be  allowed  to  vary  as  the  solution  progresses. 

For  the  case  of  subsonic  inflow,  the  stagnation  pressure,  stagnation  tempera¬ 
ture,  and  flow  angles  are  specified.  These  quantities  are  related  to  the  static  pressure 
and  static  temperature  by  the  following  equations: 


P_  _  L  _  7  ~  1  (V 


J-  /  r Tot  \ 

1  V  o.  J 


\  21 
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p_ 

pt 

u 

Vr7t 

v 

Vr7t 

w 

v™ 


_  ^  ~  ^  ( v^lY 

7  +  1  v  a.  / 

=  (  — —  )  =  specified 

\  VTot  '  0 

=  (  -  ~  )  =  specified 

V  VTot)  0 
/  w  \ 

=  1  —  )  —  specified 

V  VTot  '  0 


(6) 


The  first  two  equations  above  are  simply  the  isentropic  relations  written  in  terms 
of  the  total  velocity,  Vj-0( ,  and  the  speed  of  sound  at  a  sonic  throat,  a,.  The  speed 
of  sound  at  a  sonic  throat  is  calculated  from  the  specified  stagnation  temperature. 


K)2 


2  R  p, 
7  T  1  Pt 


Equations  6  axe  a  system  of  fb’e  equations  in  five  unknowns:  p,  p,  u,  v,  and  w, 
but  one  of  the  last,  three  equations  is  redundant.  To  complete  the  system  another 
equation  is  needed.  This  is  to  be  expected  since,  as  was  explained  earlier,  only  four 
items  may  be  specified  at  the  inflow  boundary. 

The  last,  equation  to  close  the  system  is  the  characteristic  relation  carrying 
information  upstream  to  the  inflow  boundary. 


Sp 


Su1 


t  - pc *r  “ (”'  - c) 


Sp  Su' 
6x  ~  pCJ x 


This  equation  is  forward  differenced. 


Svn 


,  (  «'  -  c)  At 

pcbuD  -  — — - Ip.'  -  Ps 


-  u'J\V 

I  ~  \  I  ’  &  /  J 


(7\ 


v  *  / 


The  subscripts  I  and  B  indicate  the  first  interior  cell  and  the  inflow  boundary  cell, 
respectively.  The  prefix  S  indicates  the  forward  in  time  difference  of  the  variable 
following  it. 

The  number  of  unknowns  is  reduced  to  three  if  the  isentropic  relations,  the 
first  two  of  Equations  6  are  written  in  terms  of  u'.  This  is  done  using  the  relation 
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(3) 


The  modified  isentropic  relations,  Equations  9,  and  the  descrete  form  of  the  up¬ 
stream  runing  characteristic  relation,  Equation  7  are  three  algebraic  relations  in 
three  unknowns. 

The  isentropic  relation  for  pressure,  the  first  of  Equation  9,  may  be  placed  in 
delta  law  form  by  considering  incremental  changes  in  the  variables  p  and  u' . 


6pg  =  Pi 


2bj 
7  +  1 


(10) 


This  equation  and  Equation  7  are  solved  for  u'.  The  pressure  and  density  are  then 
obtained  from  the  isentropic  relations,  Equations  9.  The  velocities  are  also  calu- 
lated  from  u'B  using  Equation  8  and  the  last  three  of  Equations  6.  The  specification 
of  the  flow  within  the  subsonic  inflow  boundary  cell  is  then  complete. 
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3.3.2  Subsonic  Outflow 


The  treatment  of  outflow  boundary  conditions  is  also  guided  by  the  theory  of  char¬ 
acteristics.  If  the  flow  normal  to  the  outflow  boundary  is  supersonic  then  all  char¬ 
acteristics  have  positive  slopes  and  no  information  propogates  upstream  from  the 
boundary  cell  to  the  interior  cell.  In  this  case  the  five  locally  one-dimensional  char¬ 
acteristic  equations  are  used  to  update  the  solution  in  the  boundary  cell.  If  the  flow 
normal  to  the  outflow  boundary  is  subsonic  then  the  u'-  e  characteristic  propogates 
information  upstream  from  tL°  boundary  cell  to  the  interior  cell.  In  this  case,  one 
item  must  be  specified  at  the  boundary  cell. 

For  subsonic  outflow  the  static  pressure  is  specified.  The  remaining  variables  in 
the  boundary  cell  are  calculated  using  the  four  downstream  running  characteristic 
equations  As  with  the  variables  specified  at  the  inflow  boundary,  the  requirement 
of  constant  pressure  at  the  outflow  boundary  is  written  in  delta  law  form. 


SpB  -  0 

This  equation  is  combined  with  the  four  characteristic  relations. 


(11) 


S  Oo  -j dp d 

c? 

SpB  +  pc6u'B 


Sv'B 

Sw'B 


u',At|S| 
V  oli 


PB  -  pi  +  {pB  ~  Pi) 


~  ni 


(u'  +  c)  At  |S| 
Voh 


\pb  ~  Pi  +  pc  ( u’g  -  u'j )]  =  R2 


u'rAflSI  .  ,  ,  . 

U/Aif|5  |r,  n 

-  wi\  =  R< 


The  above  four  equations  and  Equation  1 1  are  five  linear  algebraic  equations  in  the 
five  unknowns  SpB,  Su'B,  Sv'B,  Sw'B,  and  6pB.  This  system  is  solved  directly  and  the 
boundary  cell  solution  is  updated. 
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3.3.3  Noslip  Wall 


At  solid  adiabatic  walls  the  pressure  and  temperature  gradients  are  assumed  to  be 
zero  and  a  no  slip  (zero  velocity  at  the  wall)  boundary  condition  is  applied.  To 
resolve  the  boundary  layer  the  mesh  must  be  refined  near  the  wall  so  that  the  cell 
nearest  the  wall  is  deep  within  the  boundary  layer.  In  this  case,  boundary  layer 
theory  shows  that  pressure  gradient  normal  to  the  wall  is  a  higher  order  effect  and 
can  be  neglected.  Similarly,  for  adiabatic  walls  the  temperature  gradient  normal  to 
the  wall  deep  within  the  boundary  layer  is  negligable.  These  boundary  conditions 
are  satisfied  by  setting  the  pressure  and  interned  energy  in  each  boundary  cell  equal 
to  the  pressure  and  temperature  in  the  interior  cell  adjacent  to  the  boundary.  The 
no-slip  condition  is  satisfied  by  setting  the  velocity  in  the  boundary  cell  equal  and 
opposite  to  the  velocity  in  the  adjacent  interior  cell. 


3.3.4  Periodic  Boundary  Conditions 


It  is  often  possible  to  reduce  the  scale  of  a  problem  by  assuming  that  the  flow  field 
is  periodic.  For  example,  in  cascade  flows  the  cost  of  calculating  the  flow  about 
all  airfoils  of  the  cascade  may  be  prohibitive,  whereas  the  cost  of  calculating  the 
flow  about  an  isolated  airfoil  is  reasonable.  Assuming  that  the  flow  fields  about  the 
airfoils  are  identical  except  for  position,  the  cascade  problem  may  be  reduced  to  the 
calculation  of  flow  about  one  airfoil.  This  is  done  by  choosing  any  surface  passing 
above  an  airfoil  and  duplicating  it  below  the  airfoil,  as  shown  in  Figure  3,  These 
two  surfaces  define  the  solution  domain.  The  presence  of  the  neighboring  airfoils  is 
included  by  allowing  the  solution  to  vary  such  that  the  solutions  along  the  upper 
and  lower  surfaces  are  identical. 


The  periodic  boundary  conditions  are  applied  as  depicted  in  Figure  4.  On  each 
time  step  the  solution  in  the  two  layers  of  interior  cells  nearest  the  upper  boundary 
are  copied  into  boundary  cells  on  the  lower  boundary.  Likewise,  the  solution  in  the 
two  layers  of  interior  ceils  nearest  the  lower  boundary  are  copied  into  the  boundary 
cells  on  the  upper  boundary.  In  this  manner  the  solutions  extrapolated  to  tr.e  upper 
and  lower  boundaries  is  identical  and  the  solution  is  periodic. 


For  the  rotor/stator  interaction  problem  considered  in  this  investigation  the 
solution  is  periodic  in  cylindrical  coordinates.  The  flow  is  described,  however,  in 
cartesian  coordinates.  Therefore,  the  velocity  vectors  must,  be  rotated  art  "epriately 
before  they  are  placed  into  the  boundary  cells. 
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The  periodic  boundary  condition  is  very  much  like  an  interzonc  boundary  con¬ 
dition  in  that  both  involve  copying  the  solution  from  the  interior  cells  of  one  zone 
into  the  boundary  cells  of  another  zone.  For  the  periodic  boundary  conditions,  the 
zone  copied  from  and  the  zone  copied  to  just  happen  to  be  the  same  zone.  The 
implementation  of  the  periodic  boundary  conditions  actually  uses  the  same  routines 
as  the  interzonc  boundary  condition  with  the  addition  of  a  routine  to  rotate  the 
velocity  vectors.  Details  of  the  interzonc  boundary  condition  implementation  are 
given  in  Section  4. 
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I1  igure  3:  Cascade  Flow 


Figure  4:  Application  of  Periodic  Boundary  Conditions  for  cascade  flow 


4  Zonal  Boundary  Conditions 


The  pilot,  code  uses  multiple  zones  of  i,j,  k  ordered  body  grids.  There  are  two  main 
advantages  of  a  multiple  zone  code  over  a  single  zone  code.  The  first  is  that  the 
multiple  zones  offers  a  greater  degree  of  geometric  flexibility  than  is  available  with  a 
single  zone.  This  simplifies  the  task  of  generating  an  adequate  mesh  for  moderately 
complex  geometries.  For  highly  comples  geonetrics,  the  generation  of  an  adequate 
mutliple  zone  mesh  may  be  possible  when  the  generation  of  an  adequate  single  zone 
mesh  is  impossible.  The  second  advantage  is  that  the  use  of  memory  management 
is  simplified.2  The  code  may  keep  only  one  zone  within  the  primary  memory  at  a 
time  while  the  other  zones  are  in  secondary  memory  (SSD  or  disk).  For  complex 
problems  with  many  zones,  this  can  dramatically  reduce  the  amount  of  primary 
memory  required. 

The  multiple  zone  approach  has  disadvantages.  The  main  disadvantage  is  that 
flowfield  data  must  somehow  be  transferred  without  both  zones  being  present  in 
memory.  This  data  transfer  process  is  referred  to  a  specification  of  zonal  boundary 
conditions. 

The  data  transfer  at  zonal  boundaries  requires  a  twro  step  process.  The  first 
step  is  to  extract  the  flow  data  from  the  interior  of  one  zone,  the  ‘donor’  zone, 
and  store  it  in  a  temporary  location  in  memory,  called  the  ‘patch’.  At  this  point 
the  donor  zone  may  be  written  to  secondary  memory  and  the  adjacent  zone,  the 
‘host’  zone,  may  be  read  from  secondary  memory.  The  second  step  is  to  copy  the 
‘patch’  to  the  boundary  cells  of  the  ‘host,’  zone.  This  process  transfers  flow  data 
across  the  zonal  boundary  in  one  direction.  The  process  is  repeated,  with  the  zones 
trading  ‘host’  and  ‘donor’  roles,  to  transfer  flow  data  across  the  zonal  boundary  in 
the  opposite  direction. 

Two  types  of  interzone  zonal  boundary  conditions  were  utilized  during  the 
Phase  I  investigation:  a  simply  connected  butt  joint  boundary  and  a  generalized, 
sliding,  butt  joint  boundary.  The  two  differ  in  how  they  extract  data  from  the 
‘donor’  zone.  The  general  sliding  butt  joint  boundary  procedure  was  developed  as 
part  of  the  Phase  I  effort. 

2Some  form  of  memory  management  is  required  when  memory  constrained  computers  such  as 
the  CRAY  XMP  are  used 
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4.1 


Simply  Connected  Butt  Joint  Boundary 


The  procedure  for  simply  connected  butt  joint  boundaries  is  that  the  boundary 
cells  along  an  interzone  segment  of  a  zone  boundary  (i.e,,  that  segment  of  the  zone 
boundary  which  attaches  to  another  zone)  are  duplicated  in  a  patch3.  For  cases 
where  the  mesh  lines  are  continuous  across  the  zonal  interface  surface,  the  flow 
data  (conservative  field  variables,  non-conservative  field  variables,  properties,  etc.) 
are  simply  copied  from  the  appropriate  cell  of  the  donor  zone  into  the  patch.  The 
data  from  the  patch  are  then  copied  into  the  appropriate  boundary  cells  of  the  host 
zone. 


4.2  General  Sliding  Boundary 


The  primary'  goal  of  the  Phase  I  work  was  to  demonstrate  the  multi-zone  Navier- 
Stokes  analysis  on  an  unsteady  turbine  rotor/stator  flow  field.  Figure  5  shows  a 
computational  mesh  at  midspan  for  a  turbine  rotor/stator  problem.  As  shown,  the 
stator  blade  is  in  the  fixed  zone  on  the  left  side  of  the  figure  and  the  rotor  blade  is  in 
the  moving  zone  on  the  right  side.  The  two  mesh  zones  join  along  a  planar  interface. 
As  the  rotor  zone  moves  the  rotor  zone  mesh  points  slide  along  this  planar  surface. 
In  general,  mesh  points  of  the  rotor  zone  do  not  coincide  with  mesh  points  of  the 
stator  zone  at  the  “sliding  boundary” -that  is,  mesh  lines  are  discontinuous  across 
the  sliding  boundary.  This  sliding  boundary  approach  has  been  used  by  Rai[5], 
Pecry[5j,  and  others. 

The  calculation  of  fluxes  at  a  cell  face  requires  flow  field  data,  at  the  cell  centers 
of  two  cells  on  both  sides  of  the  face.  Unfortunately,  there  are  no  real  cells  for  zone 
1  that  exist  on  the  right  side  of  the  sliding  boundary  from  which  the  flux  at  cell 
faces  on  the  sliding  boundary  can  be  computed.  In  order  to  calculate  fluxes  at  the 
cell  faces  of  zone  1  on  the  sliding  boundary  flow,  field  data  from  within  zone  2  is 
interpolated  to  fictitious  boundary  cells  of  zone  1  that  overlay  zone  2  as  shown  in 
Figure  6.  Transferring  data  to  these  boundary  cells  is  difficult  for  two  reasons.  First 
of  all,  the  cell  centers  of  the  boundary  ceils  do  not  coincide  with  the  cell  centers  of 
cells  in  zone  2  and,  therefore,  require  a  three-dimensional  interpolation.  Secondly', 
in  flow  calculations  of  practical  use  each  zone  has  a  large  number  of  mesh  points. 
It  is  not  usually'  possible  to  have  both  zones  in  RAM  at  the  same  time. 

3A  separate  data  storage  location  for  the  patch  is  required  to  handle  the  situation  where  two 
large  zones  are  to  be  coupled  on  a  machine  with  limited  memory.  In  this  case,  both  zones  may  not 
fit  in  main  memory  at  the  same  time,  while  a  zone  and  a  patch  would. 
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Transferring  flow  data  across  the  sliding  boundary  was  performed  in  a  two-step 
process  as  shown  in  Figure  7.  The  first  step  consists  of  interpolating  (using  tri-linear 
interpolation)  flow  field  data  from  mesh  cells  in  zone  2  onto  the  cell  centers  of  a 
set  of  cells  that  correspond  to  the  locations  of  the  fictitious  boundary  cells  of  zone 
1.  This  set  of  cells  is  called  a  ‘patch.’  The  second  step  consists  of  simply  copying 
the  patch  data  onto  the  corresponding  boundary  condition  cells  of  zone  1.  Patch 
data  is  stored  separately  from  zone  data,  allowing  the  transfer  of  interzone  boundary 
condition  data  between  zones  with  only  one  zone  resident  in  RAM  at  a  time. 

Since  the  interpolation  must  be  repeated  every  time  step  that  the  rotor  zone 
moves,  the  interpolation  cam  become  a  major  part  of  the  computational  work.  In 
order  to  do  the  interpolation  it  is  first  necessary  to  detemine  which  cell  of  zone  2 
contains  the  center  of  each  one  of  the  patch  cells  of  zone  1.  An  exhaustive  search 
for  each  patch  cell  would  be  very  expensive.  The  search  procedure  starts  by  first 
searching  in  the  neighborhood  of  the  cell  in  zone  2  in  which  the  previous  patch  cell 
was  found.  If  this  initial  search  fails  then  the  following  seach  procedure  is  used.  The 
cells  of  zone  2  are  subdivided  into  subzones  containing  nearly  an  equal  number  of 
cells  each.  Each  of  these  subzones  is  is  in  turn  subdivided  again  into  sub-subzones. 
The  number  of  subzones  and  sub-subzones  is  variable.  The  range  coordinate  values 
in  each  coordinate  direction  of  each  subzone  (and  each  sub-subzone)  is  calculated 
at  the  beginning  of  the  time  step  prior  to  the  interpolation  process.  A  quick  search 
is  made  to  determine  in  which  subzone(s)  the  patch  cell  might  be  contained.  Next 
each  sub-subzone  of  each  selected  subzone  is  examined  to  determine  if  the  patch  cell 
lies  within  it.  Finally  after  determining  which  sub-subzone  contains  the  patch  cell 
a  cell-by-cell  search  is  performed.  This  procedure  reduces  the  time  substantially  for 
interpolation  as  compared  to  an  exhaustive  search. 

During  the  Phase  I  effort  a  form  of  trilinear  interpolation  was  used.  The  inter¬ 
polation  is  based  on  the  cell  centered  coordinates  so  that  first  step  is  to  calculate 
these  coordinates  from  the  cell  corner  point  coordinates.  For  convenience  in  the 
following  discussion,  we  call  the  point  to  which  we  are  interpolating  ‘P’.  The  next 
step  is  to  determine  which  cell  the  point  P  is  in.  The  final  step  is  to  actually  perform 
the  interpolation.  We  will  discuss  the  last  step  first  and,  to  simplify  the  discussion, 
we  will  consider  the  two-dimensional  case. 

The  interpolation  is  developed  using  a  geometric  interpretation.  As  shown  in 
Figure  8,  each  quadr  lateral  cell  in  the  mesh  may  be  divided  into  a  pair  of  triangles 
in  two  different  ways.  The  first  pair  of  triangles,  the  A-triangles,  are  obtained  by 
drawing  a  line  between  nodes  (i,j)  and  ( i  +  1  ,j  +  1).  The  second  pair  of  triangles, 
the  B-triangles,  are  obtained  by  drawing  a  line  between  the  nodes  (i,j  +  1)  and 
(i  +  1,7).  If  the  point,  P,  is  within  the  quadralateral  it  is  within  one  of  the  A- 
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Figure  5:  Midspan  Mesh  for  Turbine  Rotor/Stator  calculation 


Figure  6:  Fictitious  Boundary  Cells  of  Zone  1  Overlap  Sliding  Boundary 
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Figure  8:  Two  possible  divisions  of  a  quadralaterai  cell  into  triangles 

triangles  and  one  of  the  B-triangles.  The  first,  step  is  to  do  a  linear  interpolation 
on  each  of  the  two  triangles  that  P  is  in.  The  linear  interpolation  gives  the  value  at, 
P  as  a  weighted  average  of  the  values  at  the  vertices  of  the  triangles.  The  average 
of  the  weights  for  the  two  triangles  containing  P  give  the  weights  for  the  bilinear 
interpolation. 

The  test  to  see  if  a  point  is  within  a  triangle  is  done  by  considering  the  areas 
of  the  three  new  triangles  created  by  drawing  lines  from  P  to  the  vertices  of  the 
original  triangle.  As  shown  in  Figure  9,  P  is  outside  the  original  triangle  if  the 
sum  of  the  areas  of  the  three  new  triangles  is  greater  than  the  area  of  the  original 
triangle.  If  P  is  within  the  triangle,  the  weights  for  the  linear  interpolation  over 
the  triangle  axe  obtained  by  dividing  the  area  of  the  triangle  into  the  areas  of  the 
three  subtriangles.  The  weight  obtained  from  a  given  subtriangle  corresponds  to 
the  vertex  opposing  the  subtriangle. 

The  search  to  determine  which  quadralaterai  contains  P  is  potentially  very 
costly,  especially  when  the  total  number  of  cells  is  large.  To  simplify  the  process  the 
zone  is  broken  into  multiple  levels  of  subzones.  The  search  begins  by  determining 
which  subzone  contains  P,  then  which  sub-subzone  contains  P,  and  so  on  until 
the  subzone  which  contains  P  at  the  lowest  level  is  found.  Each  subzone  (sub¬ 
subzone,  etc.)  is  defined  by  the  smallest  horizontal  rectangle  which  envelops  the 
cells  contained  in  the  subzone,  so  the  search  through  the  subzones  is  very  efficient. 
Finally,  an  exhaustive  search  is  made  through  the  appropriate  lowest  level  subzone 
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Figure  9:  Subtriangles  with  P  within  (left)  and  outside  of  (right)  the  triangle 


until  the  actual  cell  containing  P  is  found.  The  number  of  cells  in  this  lowest 
level  subzone  (call  it  the  target  subzone)  is  relatively  small  so  the  final  search  is 
inexpensive.  The  total  cost  for  this  overall  search  is  much  less  than  the  cost  for  a 
single  exhaustive  search,  especially  when  the  number  of  mesh  points  is  large. 

The  development  of  the  trilinear  interpolation  procedure  is  similar  to  the  above 
development  of  the  bilinear'  interpolation  procedure  except  that  the  triangles  are 
replaced  by  tetrahedrons. 


5  Turbine  Rotor/Stator  Calculation 

5,1  Approach 


As  shown  in  Figure  10,  zone  1  extends  over  one  complete  period  of  the  stator  blades. 
Zone  2  extends  over  one  complete  period  of  the  rotor  blades.  For  this  analysis  the 
blade  period  of  the  rotor  blades  was  assumed  equal  to  the  blade  period  of  the 
stators.  Periodic  boundary  conditions  are  appled  between  faces  A  and  B  for  each 
zone.  Subsonic  inflow  boundary  of  zone  1.  Total  pressure,  total  temperature,  and 
flow  angle  were  specified. 


Figure  10:  Zonal  Configuration  for  Turbine  Rotor/Stator  Calculation 

As  shown  in  Figure  10  the  stator  zone  and  the  rotor  zone  are  not  circumfer¬ 
entially  aligned  at  the  sliding  boundary.  This  presents  a  problem  in  setting  flow 
field  data  in  the  interzone  patches.  Setting  the  flow  field  data  in  patch  1  consists 
of  interpolating  data  from  zone  2.  Unfortunately,  some  cells  in  patch  1  do  not  lie 
within  zone  2.  A  typical  fix  for  this  problem  would  be  to  add  some  coding  to  the 
interpolation  routine  that  would  rotate  the  patch  cell  forward  by  one  blade  period 
to  a  corresponding  position  in  zone  2  before  performing  the  interpolation.  The  in¬ 
terpolated  data  would  then  be  rotated  back  to  its  original  position.  This  method 
would  work,  but  it  would  also  complicate  the  structure  of  the  interpolation  sub- 
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routine.  Since  this  rotation  would  not  always  be  done,  a  special  test  would  have 
to  be  added  before  the  coding  that  performs  the  forward  and  reverse  rotations. 
The  interpolation  routine  would  lose  its  usefullness  for  other  applications  and  could 
introduce  bugs. 

Since  one  of  the  intents  of  this  work  was  to  develop  an  analysis  procedure 
suitable  for  general  flow  problems  involving  solid  objects  in  relative  motion,  we 
prefered  not  to  write  any  coding  that  is  specific  to  this  tost  problem.  That  is,  no 
subroutines  were  written  that  would  be  applic?  ble  to  only  the  turbine  rotor/stator 
problem.  All  additional  coding  in  the  multi-zone  3D  Navier-Stokes  analysis  must 
perform  a  simple  function,  be  of  general  use,  and  conform  to  the  layered  code 
structure  shown  in  Figure  11.  The  procedure  for  applying  interzone  boundary 
conditions  across  the  sliding  boundary  would  have  to  be  built  out  of  simpler  basic 
data  operations. 

The  present  multi-zone  3D  Navier-Stokes  analysis  is  a  result  of  an  effort  to 
develop  a  layered  code  structure.  The  planning  of  its  design  required  a  substantial 
amount  of  time.  The  development  of  the  coding,  however,  is  surprisingly  fast  be¬ 
cause  of  its  logical,  uncomplicated  code  structure.  The  layered  structure  requires 
highly  modular  programming  and  clean  data  paths  into  and  out  of  subroutines. 
This  allows  the  subroutine  modules  to  be  easily  testable  and  bugs  easily  traced  to 
offending  subroutines.  Routines  may  call  routines  in  the  same  or  any  lower  layer 
but  not  routines  in  higher  layers. 

As  shown  in  Figure  11  the  highest  level  of  layered  structure  is  the  APPLICA¬ 
TION  LAYER.  This  layer  consists  of  a  series  of  subroutine  calls.  It  is  almost  entirely 
user  definable  and  contains  some  local  processing  for  applications  specific  functions 
such  as  inputting  control  vriables.  The  subroutine  calls  from  the  APPLICATION 
LAYER  are  primarily  to  the  TASK  LAYER.  Example  TASK  LAYER  routines  are 
described  below: 

CALL  LOADZ(IZ)  Check  that  data  for  zone  IZ  is  located  in  the  program 

controlled  heap  in  RAM.  If  not,  locate  data  and  load 
it  into  RAM. 

CALL  SOLVE(IZ)  Apply  one  cycle  of  the  solution  procedure  of  the  type 

currently  assigned  to  zone  IZ,  This  would  be  one 
time-step  of  the  N.S.  solution  procedure. 

CALL  UPDATBC(IZ)  Update  all  boundary  conditions  on  zone  IZ  including 

interzone  boundary  conditions. 
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APPLICATION  LAYER 

User  Control  Program 


TASK  LAYER 

General  routines  for  initializing 
the  processing  of  data. 


LINK  LAYER 

Set  pointers  to  appropriate  data 
and  determine  data  disposition. 


ZONE  LAYER 

Processes  data.  Usually  on  a 
single  set  of  data  called  a  zone. 


UTILITY  LAYER 

Data  storage  management  (RAM  heap,  SSD,  disks), 
general  I/O,  general  data  transformation,  etc. 
Some  hardware  specific  coding. 


Figure  11:  Layered  Code  Structure  of  Pilot  Code 


TASK  LAYER  routines  primarily  call  LINK  LAYER  routines.  LINK  LAYER 
routines  look  at  various  control  data  to  determine  the  disposition  of  the  data  (in 
RAM,  SSD,  etc.)  and  the  form  of  the  data,  i.e.,  which  solver  class  (3D  N.S., 
Full-Potential,  etc.).  The  appropriate  pointer  variables  are  then  passed  on  to  the 
appropriate  routine  in  the  next  lower  layer  for  processing. 

The  routines  in  the  ZONE  LAYER  process  data.  This  includes  integration  of 
a  zone  of  field  variables  by  one  time  step,  setting  boundary  conditions  along  zone, 
boundaries  (e.g.  solid-wall  no-slip),  and  calculating  mesh  metrics  for  example.  All 
routines  in  the  ZONE  LAYER  operate  on  only  a  single  zone  or  single  set  of  data  at 
a  time. 

Data  transfer  between  zones  is  handled  as  described  previously  in  a  two-step 
procedure  using  ‘patches.’  ZONE  LAYER  routines  exist  for:  1)  interpolating  data 
from  one  zone  into  cells  of  a  patch,  and  2)  copying  patch  data  into  zone  boundary 
condition  data.  The  patch  serves  as  a  buffer  for  temporary  data  storage.  The  ZONE 
LAYER  routines  can  be  relatively  simple  to  write,  since  each  routine  operates  on 
oniy  one  zone  at  a  time.  The  problems  in  managing  the  multi-zone  data  are  handled 
by  routines  in  the  APPLIATION,  TASK,  and  LINK  LAYERS. 

The  lowest  layer  is  called  the  UTILITY  LAYER..  The  routines  in  this  layer 
manage  the  computer  resources.  The  perform  operations  that  are  sometimes  ma¬ 
chine  specific.  These  include  1)  initializing  the  RAM  heap,  requesting  and  returning 
blocks  of  data  to  the  heap,  2)  moving  zone  and  patch  data  between  RAM,  SSD, 
and  disks,  and  3)  generalized  input  and  output  of  user  data. 

The  procedure  for  applying  the  interzone  boundary  condition  is  being  developed 
with  the  layered  structure  in  mind.  The  layered  structure  separates  the  complexities 
of  managing  multiple  zones  of  data  from  the  nitty  gritty  low-level  date  manipula¬ 
tions.  This  simplifies  the  coding  greatly,  but  must  be  well  thought  out  to  have 
general  application.  The  present  layered  structure  allows  the  interne  boundary 
condition  to  be  applied  very  easily.  To  do  this,  the  two  following  basic  routines  are 
added  to  the  TASK  LAYER,  (plus  the  auxiliary  lower-level  routines  that  perform 
the  data,  manipulation). 

DUPZON(IZ,IDZ)  This  routine  creates  a  duplicate  copy  of  zone  IZ  and 

assignes  it  the  zone  number  IDZ. 

ROTZON(IZ,ANG).  This  routine  rotates  the  coordinates  and  vector  com¬ 

ponents  about  an  axis  passing  through  the  origin  in 
the  Cartesian  x-direction. 
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Figure  12:  lnterzone  Data  Transfer  Between  Non-aligned  Zones 

Setting  the  interzone  boundary  conditions  across  the  sliding  boundary  becomes 
trivial  with  the  use  of  the  above  routines.  Consider  setting  the  interzone  boundary 
conditions  on  zone  1.  The  first  step  is  to  interpolate  flow  field  data  from  zone  2  into 
patch  1  (which  is  associated  with  the  boundary  cells  of  zone  1).  Prior  to  setting  the  : 

cel's  in  patch  1,  a  third  zone  is  temporarily  created  by  duplicating  zone  2  using  the  j 

TASK  LAYER  routine  DUPZON.  This  temporary  zone  is  given  the  number  3  and  ( 

rotated  counterclockwise  one  blade  period  so  that  it  trails  behind  zone  2  as  shown 
in  Figure  12.  J 

Setting  data  into  patch  1  is  now  straight  forward  and  proceeds  as  follows.  Zone 
2  and  patch  1  are  loaded  into  RAM.  Flow  field  data  is  interpolated  from  zone  2  into 
the  ceils  of  patch  1  that  lie  within  zone  2.  Next,  zone  2  is  moved  from  RAM  to  the 
SSL*  and  zone  3  is  loaded  into  RAM.  Flow  field  data  in  zone  3  is  interpolated  into 
the  remaining  cells  in  patch  1  which  lie  within  zone  3.  Zone  3  is  then  deleted.  As 
shown  in  Figure  12,  the  process  occurs  in  a  similar  manner  for  setting  cells  in  patch  t 

2  for  she  boundary  conditions  of  zone  2. 

The  application  of  periodic  boundary  conditions  on  zones  1  and  2  was  also 
performed  with  the  use  if  simple  TASK  LAYER  routines.  Zone  1  is  shown  in 
isolation  in  Figure  13.  Data  from  two  planes  of  internal  cells  parallel  to  face  A 
are  transferred  into  patch  4  using  a  TASK  Layer  routine  called  UPDZIP(IP,IZ). 
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Figure  13:  Periodic  Boundary  Conditions  for  Stator  Zone 

TJPDZIP(IP,IZ)  extracts  data  from  zone  1  (IZ=1)  and  copies  it  into  patch  4  (IP=4). 
Another  TASK  LAYER  routine  called  ROTPCH(IP,IZ,ANG)  then  rotates  patch  4 
(IP=4)  by  an  angle  of  ANG  about  an  axis  through  the  origin  in  the  Cartesian  x- 
direction.  Setting  ANG  equal  to  the  negative  of  the  angular  blade  period  will  rotate 
the  patch  (and  vector  components  in  its  data)  so  that  it  coincides  with  the  boundary 
condition  cells  labeled  B  in  Figure  13.  Then  using  a  third  TASK  LAYER  routine 
called  UPDIPZ(IP,IZ)  data  transferred  from  patch  4  (IP— 4)  to  the  boundary  cells 
in  zone  (IZ=1). 

The  above  TASK  LAYER  routines  DUPZON  and  ROTZON  can  be  used  in 
other  applications  as  well.  They  are  now  building  blocks  within  a  library  of  TASK 
LAYER  routines.  They  do  not  have  to  be  bugged  again.  Once  a  large  library  of 
TASK  LAYER  routines  have  been  built,  the  user  can  create  applications  to  many 
flow  problems  without  coding  a  single  low-level  routine.  This  layered  approach 
also  reduced  the  probability  of  programming  errors,  since  the  low-level  routines  are 
relatively  simple  functions  that  can  be  easily  tested,  the  complexities  of  data  and 
control  management  are  confined  to  the  high-level  layers,  and  new  applications  may 
require  only  changes  to  routines  at  the  TASK  LAYER. 

A  significant  contribution  of  the  Phase  I  work  is  the  demonstration  of  the  ben¬ 
efits  of  a  layered  code  structure.  The  development  of  a  sophisicated  well-engineered 
protocol  for  a  layered  code  structure  would  be  of  great  usefulness  to  CFD  users 
and  developers  of  flow  analysis  programs  and  other  physical  process  simulation 
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programs. 


This  layered  approach  has  a  number  of  benefits:  1)  it  greatly  reduces  probabil¬ 
ity  of  programming  errors,  2)  substantially  reduces  programmimg  time  in  applying 
codes  to  new  problems,  3)  has  a  potential  for  becoming  a  programming  standard 
allowing  networked  users  to  benefit  from  other  users’  programming  (i.e.  share  rou¬ 
tines),  and  4)  allows  researchers  interested  in  algorithm  development  a  convenient 
environment  to  conduct  research  wothout  spending  time  developing  a  complete  pro¬ 
gram. 


5.2  Results  of  Calculations 


The  above  code  is  being  applied  to  the  solution  of  a  turbine  rotor/statoi  interaction 
problem.  Geometry  is  based  on  the  rctor/stator  model  tested  experimentally  by 
Dring  et  al.  [4],  Calculations  of  this  flow  have  beer  presented  by  Rai[l,3].  Toe 
model  tested  experimentally  by  Dring  had  22  stator  blades.  Calculations  of  this 
flow  have  been  presented  by  Rai[l,3].  The  model  tested  experimentally  by  Dring 
had  22  stator  blades  and  28  rotor  blades.  A  calculation  using  the  model  geometry 
would  require  that  11  stator  blades  and  14  rotor  blades  be  modeled  before  periodic 
boundary  conditions  could  be  used.  To  mimimizc  the  computational  expense  we 
make  the  same  assumption  as  Rai[3]  and  model  the  case  with  22  rotor  and  22  stator 
blades  by  enlarging  the  rotor  blade.  No  attempt  is  being  made  to  model  the  rotor 
blade  tip  effects. 

Two  calculations  were  performed.  The  first  was  with  the  rotor  blade  stationary 
and  the  second  was  with  the  rotor  blade  moving.  In  both  calculations  the  mesh 
lines  were  discontinuous  across  the  boundary  between  the  rotor  and  stator  zones, 
and  the  general  sliding  zonal  boundary  condition  was  used.  The  first  calculation, 
with  the  rotor  fixed,  provided  an  excellent  test  of  the  zonal  boundary  condition, 
without  the  added  complexity  of  the  moving  mesh. 


The  mesh  used  for  the  stationary  rotor  calculation  is  shown  in  Figure  14.  An 
O-mesh  is  used  for  both  the  stator  and  the  rotor.  The  dimensions  of  the  mesh  are 
66x23x4  for  the  stator  zone  and  74x23x4  for  the  rotor  zone.  The  initial  conditions 
for  the  calculation  were  a  uniform  axial  flow  of  100  m/s.  The  flow  at  the  inflow 
boundary  was  specified  as  axial  with  a  stagnation  pressure  of  1.0  MPa  and  a  stag¬ 
nation  temperature  of  1000.0  degrees  Kelvin.  The  outflow  pressure  was  specified 
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Figure  14:  Mesh  used  for  stationary  rotor  calculation 

as  0.896  MPa.  The  solution  was  obtained  after  8100  explicit  time  steps  4.  At  this 
point  the  L2-norm  of  the  residual  had  dropped  four  orders  of  magnitude.  The  pres¬ 
sure  contours  at  the  median  radial  surface  (midspan)  are  shown  in  Figure  15,  the 
velocitiy  vectors  are  shown  in  Figure  16,  and  the  streamlines  are  shown  in  Figure 
17.  These  results  are  qualitatively  correct. 

The  mesh  used  for  the  moving  rotor  calculation  is  shown  in  Figure  18  at  three 
times  dining  the  cycle.  A  coarse  mesh  was  used  for  this  calculation  because  of 
limited  computer  resources.  Shortly  after  the  last  time  in  Figure  18  the  rotor  zone 
is  rotated  back  5  and  the  solution  is  continued.  The  velocity  vectors  at  these  three 
times  are  shown  in  Figure  19  and  the  streamlines  at  the  first  two  times  are  shown 
m  Figure  20. 


4 Local  time  stepping  was  used  to  accelerate  the  convergence. 

5The  physical  interpretation  is  that  the  calculation  proceeds  with  the  next  rotor  in  line. 
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6  Conclusions 


The  goal  of  the  Phase  I  effort  was  to  demonstrate  tiie  feasibility  of  calculating  flows 
with  bodies  in  relative  motion.  With  this  goal  in  mind,  the  results  of  the  Phase  I 
effort  are  encouraging.  It  has  shown  the  benefits  of  developing  a  layered  software 
structure,  and  demonstrated  the  feasibility  of  the  coupling  procedure  for  interzone 
boundary  conditions  between  zones  in  relative  motion.  It  is  felt  that  the  procedures 
developed  during  Phase  I  can  form  the  basis  of  a  versatile  code  for  calculating 
complex  flows  containing  bodies  in  relative  motion. 
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7  Recommendations 


Since  the  primary  objectives  of  the  Phase  I  work  were  achk  ,i,  the  following  rec¬ 
ommendations  are  made  for  follow  up  work.  The  ultimate  goal  of  this  effort  is 
to  develop  a  versatile  computer  program  for  calculating  complex  flows  containing 
bodies  in  relative  motion. 


1.  Software  Architecture:  Further  study  is  needed  to  determine  the  optimal  soft¬ 
ware  architecture.  The  objective  is  to  develop  a  software  architectural  plan  for 
building  an  integrated  modular  analysis  system  that  is  capable  of  efficiently 
modeling  the  complex  physical  interactions  which  occur  when  neighboring 
bodies  are  in  relative  motion.  The  data  structure  of  the  analysis  system  will 
be  zonal  in  nature.  Each  zone  may  have  a  different  solver  class  and,  for  that 
matter,  a  completely  different  data  structure.  For  instance,  the  zones  used  to 
model  the  fluid  dynamics  (Navier- Stokes,  Euler,  full-potential  equations)  may 
have  a  structured  grid  (a  grid  with  a  logical  i,j,k  ordering),  and  the  zones  used 
to  model  structural  dynamics  may  have  an  unstructured  finite-element  grid. 

2.  Multiple  Solvers:  There  are  many  advantages  to  incorporating  many  dissim¬ 
ilar  solvers  into  a  single  code.  The  engineer  needs  to  learn  only  one  user 
interface  and  the  various  physical  models  can  share  resources  such  as  plotting 
routines,  I/O  routines,  linear  algebra  routines,  etc.  In  addition,  it  facilitates 
the  coupling  of  different  physical  models,  such  as  a  fluid  dynamics  and  struc¬ 
tural  dynamic  solvers.  The  successful  integration  of  widely  differing  solution 
procedures  into  a  single  program  requires  careful  consideration  of  the  soft¬ 
ware  architecture,  as  described  in  Section  5.  A  layered  software  architecture, 
such  as  that  developed  during  Phase  I,  is  a  likely  candidate  for  isolating  the 
difference  between  solvers. 

3.  Zonal  Coupling:  The  next  step  is  to  couple  the  different  solvers  described 
above.  This  would  allow  complex  physical  interactions,  such  as  fluid/structural 
interaction,  to  be  studied.  Again,  the  goal  would  like  to  keep  the  program  as 
versatile  as  possible,  so  that  it  would  not  be  limited  to  any  one  category  of 
problems. 


Amtec  Engineering,  Inc.  feels  that  the  development  of  such  versatile  analysis  soft¬ 
ware  would  benefit  the  U.S.  Army  and  the  engineering  community  as  a  whole. 


40 


References 


[lj  Anderson,  D.A.,  Tannehill,  J.C.,  and  Fletcher,  R.H.  Computational  Fluid.  Me¬ 
chanics  and  Heat  Transfer ,  McGraw-Hill,  1984. 

[2]  Baldwin,  B.  and  Lomax,  H.,  “Thin-layer  Approximation  and  Algebraic  Model 
for  Separated  Turbulent  Flows,”  AIAA  PAper  No.  78-257,  Jan.  1978. 

[3]  Steger,  J.L.  and  Warming,  R.F.,  “Flux  Vector  Splitting  of  the  Inviscid  Gasdy- 
namic  Equations  with  Applications  to  Finite-Difference  Methods,”  Journal  of 
Comp.  Phys.,  Vol.  40,  pp263-293,  1981. 

[4]  Pulliam,  T.H.  and  Seger,  J.L.,  “On  Implicit  Finite  Difference  Simulations  of 
Three-Dimensional  Flows”  AIAA  Paper  No.  78-10,  Jan.  1978. 

[5]  Rai,  M.M.  “Navier-Stokes  Simulations  of  Rotor-Stator  Interaction  Patched  and 
Overlaid  Grids,”  AIAA  Paper  No.  85-1519,  July  1985. 

[6]  Peery,  K.M.,  Ponten,  B.D.,  and  Roberts,  D.W.,  “Simulation  of  Unsteady  Two- 
Dimensional  Inviscid  Flow  Fields  Around  Geometrically  Complex  Objects,” 
AIAA  Paper  No.  85-1273. 

[7]  “Unsteady  Three  Dimensional  Navier-Stokes  Simulations  of  Turbine  Rotor- 
Stator  Interactions  Including  Tip  Effects,”  AIAA/ASME/SAE  23rd  Joint 
Propulsion  Conference,  San  Diego,  CA,  June  1987. 

[8]  Dring,  R.P.,  Joslyn,  H.D.,  Hardin,  L.W.,  and  Wagner,  J.H.,  “Turbine  Rotor 
Stator  Interaction,”  J.  of  Engrg.  for  Power ,  Vol.  104,  Oct.  1982. 

[9]  Chima,  R.V.,  “Development  of  an  Explicit  Multigrid  Algorithm  for  Quasi- 
Three-Dimensionai  Viscous  flows  in  Turbonmaclunery,”  NASA  TM-87128,  186. 

[10]  Davis,  R..L.,  Ni,  R.H.,  and  Carter,  J.E.,  “Cascade  Viscous  Flow  Analysis  Using 
the  Navier-Stokes  Equations,”  AIAA  Paper  No.  86-0033,  Jan.  1986. 

[11]  Pulliura,  T.H.  and  Steger,  J.L.  “On  Implicit  Finite  Difference  Simulations  of 
Three  Dimensional  Flows,”  AIAA  Paper  No.  78-10,  January  1978. 


